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In  this  work,  a  three-dimensional,  non-isothermal  and  two-phase  computational  fluid  dynamics  model  of 
a  proton  exchange  membrane  (PEM)  fuel  cell  with  straight  flow  field  channel  is  developed  and  validated. 
The  model  is  used  to  predict  the  performance  of  the  PEM  fuel  cell  with  changing  parameters  of  the  cathode 
catalyst  layer  which  was  usually  assumed  to  be  composed  of  spherical  agglomerates.  The  effect  of  cathode 
catalyst  layer  parameters  such  as  catalyst  layer  thickness,  ionomer  film  thickness,  agglomerate  size  and 
porosity,  on  the  current  density  and  power  output  of  the  PEM  fuel  cell  is  investigated.  The  numerical 
results  reveal  that  competitive  influence  of  resistances  to  transport  of  species,  electron  and  proton  within 
the  cathode  catalyst  layer  determines  the  performance  of  the  PEM  fuel  cell  in  terms  of  area  specific  power 
density  (W  cm-2 )  and  mass  specific  power  density  (kW  g"1 ). 

©  2010  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

In  the  last  decade,  proton  exchange  membrane  (PEM)  fuel  cells 
have  received  considerable  attention  due  to  their  high  efficiency, 
low  temperature  operation  range,  scalability,  compactness  and 
low  emissions  [1,2].  However,  there  are  many  technical  barriers 
that  prevent  complete  commercialization  of  PEM  fuel  cells.  Among 
these  reduction  of  precious  metal  catalyst  loading,  reliability  and 
stability  of  performance  under  variety  of  conditions  and  durabil¬ 
ity  of  fuel  cell  module  components  are  the  most  important  issues 
that  significant  research  effort  has  been  devoted  [3-5].  Catalyst 
loading  and  reliability,  stability  and  durability  of  PEM  fuel  cell 
modules  may  be  improved  by  optimizing  the  design  and  operat¬ 
ing  conditions.  This  requires  complete  understanding  of  complex 
and  coupled  transport  phenomena  of  mass,  energy  and  current. 

It  is  well  known  that  the  performance  of  the  PEM  fuel  cells 
is  strongly  dependent  on  the  electrodes,  especially  at  the  cath¬ 
ode.  The  critical  component  of  the  electrodes  is  the  catalyst  layers. 
The  catalyst  layers  are  typically  composed  of  carbon,  ionomer  and 
platinum  which  are  randomly  dispersed  within  a  porous  matrix. 
This  random  structure  makes  it  difficult  to  find  the  optimum  cata¬ 
lyst  layer  composition  as  well  as  structure.  However,  a  properly 
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constructed  mathematical  model  with  advanced  manufacturing 
techniques  might  help  engineers  to  obtain  the  highest  possible 
performance  with  minimum  use  of  precious  metal  catalyst. 

In  the  last  decade,  the  agglomerate  structure  model  for  cata¬ 
lyst  layers  has  received  a  lot  of  attention  in  modeling  of  the  PEM 
fuel  cells  [6-12].  This  is  because  of  various  experimental  studies 
supporting  the  presence  of  agglomerates  and  modeling  studies  in 
which  the  prediction  capability  of  the  agglomerate  model  revealed 
[6,10-13].  Agglomerate  model  assumes  that  platinum  supported 
carbon  particles  are  grouped  and  bonded  with  ionomer  to  form 
the  agglomerates.  These  agglomerates  may  also  be  surrounded 
by  a  thin  ionomer  film.  The  shape  of  the  agglomerates  is  usually 
assumed  to  be  spherical  since  the  spherical  agglomerate  model  is 
showed  to  be  the  realistic  representative  of  catalyst  layers  [10]. 
For  the  spherical  agglomerate  model,  the  electrochemical  reactions 
are  usually  modeled  with  well  known  internal  effectiveness  factor 
approach  commonly  used  for  heterogeneous  reactions  occurred  in 
porous  catalyst  particles. 

In  the  literature,  various  numerical  parametric  studies  were 
presented  for  analyzing  the  effect  of  the  cathode  catalyst  layer 
parameters  on  the  PEM  fuel  cell  performance  [14-20].  To  the  best 
knowledge  of  the  authors  of  this  article,  none  of  these  studies  are 
performed  for  non-isothermal,  two-phase  operation  of  a  PEM  fuel 
cell  using  three-dimensional  geometry  and  spherical  agglomer¬ 
ate  model  for  cathode  catalyst  layer.  However,  in  the  PEM  fuel 
cells  three-dimensional  transport  effects  on  distribution  of  species 
concentration,  temperature,  etc.,  is  not  negligible.  Hence,  in  this 
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Nomenclature 

Ea 

activation  energy  (J  mol-1 ) 

aagg 

active  agglomerate  surface  area  per  unit  volume  of 
agglomerate  (m-1) 

teat 

catalyst  layer  thickness  (fxm) 

mpt 

catalyst  loading  per  unit  MEA  area  (mgPt  cm-2) 

C 

concentration  (kmol  m-3) 

kc 

condensation  rate  coefficient  (s-1 ) 

J 

current  density  (A cm-2)  or  (A  mg"1 ) 

P 

density  (kgm-3) 

D 

diffusion  coefficient  (m2  s-1 ) 

x1 

1 

QJ 

b 

diffusion  layer  electric  conductivity  in  x-  and  y- 
direction  (Snrr1) 

G e-z 

diffusion  layer  electric  conductivity  in  z-direction 
(Sm-1) 

rPt 

effective  catalyst  surface  area  ratio 

Peff 

effectiveness  factor 

O 

electric  or  proton  conductivity  (S  m-1 ) 

nd 

electro-osmotic  drag  coefficient 

fee 

evaporation  rate  coefficient  (Pa-1  s-1 ) 

F 

Faraday  constant  (C  mol-1 ) 

hc 

heat  of  condensation  of  water  (J  kg-1 ) 

Ho2 

Henry’s  law  constant  for  the  dissolution  of  oxygen 
into  the  electrolyte  (Pa  m3  kmol"1 ) 

S 

ionomer  film  thickness  (nm) 

I< 

isotropic  permeability  (m2) 

k 

isotropic  thermal  conductivity  (W  m-1 I<-1 ) 

EWm 

membrane  equivalent  weight  (kg  kmol"1 ) 

M 

molecular  weight  (kg  kmol"1 ) 

V 

over-potential  (V) 

Po2 

partial  pressure  of  oxygen  (Pa) 

I<xy 

permeability  inx-  andy-direction  (m2) 

I<z 

permeability  in  z-direction  (m2) 

rPtC 

platinum  mass  ratio  on  Pt|C  catalyst 

£0 

porosity  (volume  fraction  of  voids)  of  fuel  cell  com¬ 
ponents 

A  HZ 

reaction  enthalpy  for  anode  half  reaction 
(J  kmol-1  K-1 ) 

AH- 

reaction  enthalpy  for  cathode  half  reaction 
(Jkmol-1  K-1) 

Jo 

reference  exchange  current  density  (Am-2) 

koRR 

reaction  rate  coefficient  for  oxygen  reduction  reac¬ 
tion  (ms-1) 

s 

source 

s 

saturation 

Yd 

saturation  exponent  for  diffusion  coefficient  correc¬ 
tion 

k-dw 

sorption/desorption  rate  coefficient  (s-1) 

Sc 

specific  surface  area  of  catalyst  particles  (m2  g-1 ) 

£agg 

spherical  agglomerate  porosity  (volume  fraction  of 
voids  in  the  spherical  agglomerates) 

ragg 

spherical  agglomerate  radius  (p>m) 

kxy 

thermal  conductivity  in  x-  and  y-direction 
(Wm-1 1<-1) 

kz 

thermal  conductivity  in  z-direction  (W  m-1  K-1 ) 

V 

Thiele  modulus 

a 

transfer  coefficient 

u 

velocity  x  component  (m  s-1 ) 

V 

velocity y  component  (ms-1) 

w 

velocity  z  component  (m  s-1 ) 

p 

viscosity  (kg  m-1  s-1 ) 

£i-film 

volume  fraction  ionomer  covering  the  spherical 
agglomerates 

Si  volume  fraction  ionomer  in  the  catalyst  layer 

(£i-agg  +  £ i-film ) 

£i-agg  volume  fraction  ionomer  in  the  spherical  agglomer¬ 
ates 

sS0i  volume  fraction  of  solid  phase 
Sft  volumetric  transfer  current  (A  m-3 ) 

X  water  content  (kmol  H2O  kmol-1  SO3) 

Superscripts 
CL  catalyst  layer 

DL  diffusion  layer 

eff  effective  property 

MP  mono-polar  plate 

ref  reference  value 

Subscripts 
a  anode 

C  carbon 

c  cathode 

dw  dissolved  water 

e  electric 

i  ionomer 

M  membrane 

Pt  platinum 

p  proton 

sol  solid 


study,  a  three-dimensional,  non-isothermal  and  two-phase  compu¬ 
tational  fluid  dynamics  model  is  used  for  systematic  analysis  of  the 
performance  of  a  PEM  fuel  cell  by  considering  area  specific  and  plat¬ 
inum  mass  specific  polarization  curves.  In  the  model,  ionomer  film 
covered  spherical  agglomerate  structure  is  used  for  the  modeling 
of  cathode  catalyst  layer. 

2.  Three-dimensional  PEM  fuel  cell  model 

A  typical  PEM  fuel  cell  consists  of  a  five-layer  membrane- 
electrode  assembly  (MEA)  sandwiched  in  between  two  mono-polar 
plates  (or  two  bi-polar  plates  in  a  stack)  on  which  anode 
and  cathode  flow  fields  are  grooved.  A  five-layer  MEA  con¬ 
tains  two  electrodes  as  anode  and  cathode  separated  by  a 
proton  conducting  polymer  electrolyte  membrane  (PEM).  Each 
of  these  electrodes  consists  of  a  thin  catalyst  layer  coated  on 
one  side  of  diffusion  layer  which  serves  as  diffusion  medium  for 
transport  of  reactant  and  products  as  well  as  for  transport  of 
electrons. 

In  a  typical  PEM  fuel  cell  operation,  usually  humidified  hydro¬ 
gen  gas  is  supplied  to  anode  flow  channel  which  diffuses  through 
anode  diffusion  layer  until  reaching  the  anode  catalyst  layer  where 
protons  (H+)  and  electrons  are  formed  at  the  catalyst  surface. 
The  protons  are  transferred  through  ionomer  inside  the  cata¬ 
lyst  layer  and  the  membrane  to  the  cathode  catalyst  layer  and 
electrons  are  transferred  through  the  external  circuit  to  cathode 
catalyst  layer.  Meanwhile,  humidified  oxygen  (or  air)  is  fed  to 
the  cathode  flow  channel  and  oxygen  reaches  the  cathode  cat¬ 
alyst  layer  by  diffusing  through  the  cathode  electrode.  At  the 
catalyst  surface,  protons,  electrons  and  oxygen  forms  water.  The 
net  reaction  in  the  PEM  fuel  cell  is  given  by  following  reac¬ 
tion: 

H2  +  \  O2  H20  +  Heat  +  Electricity  (1) 
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Fig.  1.  Schematic  representation  of  volume  fraction  of  each  phase  in  the  PEM  fuel 
cell  catalyst  layer. 


2.2.  Model  assumptions 

The  following  key  assumptions  are  made:  (a)  steady-state,  non- 
isothermal  operation  of  PEM  fuel  cell  is  considered;  (b)  catalyst 
layers  are  assumed  to  be  isotropic  and  homogeneous;  (c)  identical¬ 
sized,  evenly  dispersed  spherical  agglomerates  are  used  to  describe 
catalyst  layers;  (d)  each  spherical  agglomerate  is  covered  by  a  thin 
constant  thickness  ionomer  film;  (e)  it  is  assumed  that  voids  in 
the  spherical  agglomerates  filled  only  with  ionomer;  (f)  gases  are 
assumed  to  be  ideal,  and  gas  flow  is  laminar;  (g)  water  is  produced 
in  the  ionomer  phase  of  the  cathode  catalyst  layer;  (h)  anisotropic 
properties,  such  as  thermal  conductivity,  is  considered  for  the  dif¬ 
fusion  layers;  (i)  fine  mist  flow  is  assumed  for  liquid  water  in  gas 
channels. 


2.2.  Structure  of  the  cathode  catalyst  layer 


In  the  cathode  catalyst  layer  solid,  gas,  liquid  and  ionomer  vol¬ 
ume  fractions  must  be  properly  defined.  In  Fig.  1  volume  fraction  of 
each  phase  is  shown  schematically  on  unit  catalyst  layer  volume. 

The  value  of  solid  phase  volume  fraction  (. £sol )  can  be  calculated 
using  platinum  loading  (mPt),  catalyst  layer  thickness  (tcat)  and  Pt|C 
ratio  (rPtc): 

-  -  (s  ■ +  )  £ 

Ionomer  volume  fractions  are  given  by 


i-film 


£ sol 


(1  —  £agg ) 


(3) 


'SO 


l£agg 


'i-agg 


' agg 


(4) 


where  £agg  is  the  porosity  or  the  volume  fraction  of  voids  between 
solid  particles  of  the  spherical  agglomerates.  The  total  ionomer 
volume  fraction  of  the  catalyst  layer  £,•  is  simply  the  summation 
of  Si-fiim  and  £i-agg •  For  two-phase  operation  gas  phase  volume 
fraction  is  obtained  using  liquid  water  saturation  (s)  and  the  dry 
porosity  (s0 )  that  is  the  volume  fraction  of  voids  in  the  porous  layers 


£gas  —  £()(1  —  5)  (5) 

The  active  surface  area  of  catalyst  particles  per  unit  agglomerate 
volume  ( aagg )  has  to  be  known  to  calculate  the  volumetric  transfer 


current  in  terms  of  catalyst  volume.  Thus  following  expression  is 
used 


dagg 


=  Scr  Pt 


mPt 

feat  (1 


1 

£0  ~  £i-film ) 


(6) 


where  Sc  is  the  specific  surface  area  of  the  catalyst  particles.  Fol¬ 
lowing  correlation  obtained  from  the  data  reported  by  the  catalyst 
manufacturer  E-TEK  is  used  in  the  present  model  [21  ] 

Sc  =  7.401  X  106rir  -  1.811  X  lO7^ r 


+  1.545  x  107rj;  -6.453  x  106rPtc  +  2.054  x  106  (7) 


Specific  area  of  the  catalyst  particles  is  multiplied  by  an  effective 
surface  area  ratio  (rPt)  in  order  to  take  into  account  the  isolation  of 
active  area  [6,22]. 

2.3.  Governing  equations  and  source  terms 

The  governing  equations  of  the  present  model  which  are  listed  in 
Table  1  include  conservation  of  mass,  momentum,  gas  species  (H2, 
02  and  H20),  energy,  electric  and  proton  potentials,  dissolved  water 
and  liquid  water.  All  governing  equations  are  listed  in  conservative 
form,  appropriate  for  use  in  commercial  flow  solvers. 

In  the  porous  zones  of  the  PEM  fuel  cell,  effects  of  porosity  and 
liquid  water  are  taken  into  account  in  the  source  terms  of  Eqs. 
(9)— ( 11).  Effective  diffusion  coefficients  (Dy^)  of  individual  species 
are  calculated  accordingly  to  include  the  effect  of  porous  media, 
through  Bruggemann’s  correction  [6,12,23],  and  liquid  water  [24]: 

of  =  e™0-s)nDj  (20) 

where  eg,  s,  yd  and  Dj  are  the  gas  phase  volume  fraction,  liquid 
water  saturation,  saturation  dependence  exponent  and  diffusion 
coefficient  of  gas  component  j  in  gas  mixture,  respectively.  In  the 
gas  channels,  effective  diffusion  coefficient  reduces  to  diffusion 
coefficient  of  component  j  in  the  gas  mixture  (Dj). 

Consumption  and  production  of  gas  species  are  defined  in  the 
source  terms  of  species  conservation  equations  (Eqs.  (12)-(14)). 
Mass  exchange  of  water  between  gas  and  liquid  phases  is  included 
in  the  source  term  of  Eq.  (14)  for  catalyst  and  diffusion  layers.  In 
the  catalyst  layers  source  term  also  includes  sorption/desorption 
of  water.  Conservation  of  energy  is  represented  by  Eq.  (15)  where 
heat  source  term  {Sth )  includes  electrochemical  reaction  heat  ( Srxn ), 
ohmic  heat  (Sohm_e  and  Sohm_i)  and  phase  change  heat  terms; 
condensation/evaporation  (S/at)  and  sorption/desorption  (Ssd)  heat 
sources.  The  heat  of  sorption/desorption  was  assumed  to  be  equal 


Table  1 

Governing  model  equations. 


Conservation  equation 

Conservative  form 

Mass 

^  ■  (Pgh)  —  Smass 

(8) 

Momentum—  x 

V  •  ( PgVU )  =  -dP/dx  +  V  •  (/Xg  Vli)  +  Smom-x 

(9) 

Momentum— y 

V  ■  ( pgVV )  =  -dP/dy  +  V  ■  (/XgVy)  +  Smom_y 

(10) 

Momentum—  z 

V  ■  (PgUW)  =  - 3P / 3 Z  +  V  •  (/Xg  Vw)  +  Smom-z 

(11) 

Species-Eh 

V  ■  (pgvYH2 )  =  V  •  (D|  VYh2)  +  Vh2 

(12) 

Species-02 

V  •  (pgvYo2 )  =  v  ■  (Df2  V Yo2 )  +  Ssp.o2 

(13) 

Species-E^O 

V  ■  (PgUYH2o)  =  V  •  (^h^o^h2°)  +  Ssp.H20 

(14) 

Energy 

V  •  ( pgCpvT )  =  v  •  (fc^vr)  +  Sth 

(15) 

Electric  potential 

0  =  V(CTfV0e)  +  Se 

(16) 

Proton  potential 

0  =  V -(of  V<pp)  +  Sp 

(17) 

Dissolved  water 

0  =  V  •  (Ddw V Cdw)  +  £dw 

(18) 

Liquid  water 

V  •  (pif(s)v)  =  V  ■  (ptDsVs)  +  Ss 

(19) 
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to  the  latent  heat  of  water  evaporation/condensation  [27].  Trans¬ 
ports  of  electron  and  proton  were  modeled  as  in  Eqs.  (16)  and  ( 1 7)  in 
which  effective  electric  and  proton  conductivities  were  calculated 
using  Bruggemann’s  type  correction: 

CFf  =  s'J(re  (21) 

of  =  e]'5<yP  (22) 

In  the  present  model,  anisotropic  electrical  conductivity  values 
are  used  for  diffusion  layers.  Proton  conductivity  (crp)  of  ionomer 
phase  in  catalyst  layers  and  the  membrane  were  calculated  using 
the  expression  given  in  the  work  of  Springer  et  al.  [25] 


gp  =  (0.51397. 


0.326)  exp 


1268 


303 


X  >  1 


(23) 


The  conservation  equation  of  dissolved  water  in  the  ionomer 
phase  contains  diffusion  term  and  electro-osmotic  drag  terms, 
latter  was  included  as  a  source  term  in  Eq.  (18).  The  diffusion 
coefficient  of  dissolved  water  in  membrane  is  calculated  using  the 
correlation  given  by  Motupally  et  al.  [26] 

3.1  x  10_7A(exp(0.28A)  -  l)exp(-2346/7)  0  <  A  <  3 

4.17  x  10“8A(1  +  161  exp(-A))exp  (-2346/r)  3  <  A  <  17 

Liquid  water  governing  equation  is  given  in  Eq.  (19)  where 
viscous  drag,  capillary  diffusion  and  phase  change  terms  were 
included  [22].  In  the  porous  diffusion  and  catalyst  layers,  capillary 
pressure  induced  diffusion  term  dominates  since  viscous  drag  term 
is  insignificant  due  to  small  gas  velocities  inside  the  porous  media. 
Liquid  water  diffusivity,  also  called  as  capillary  diffusivity,  in  Eq. 
(19)  was  calculated  using  [23]: 


Kj  dpc 
fii  ds 


(25) 


where  I<t  is  the  liquid-phase  permeability  and  pc  is  the  capillary 
pressure.  Both  are  functions  of  liquid  water  saturation  (s)  and  are 
defined  as 


Ki(s)  =  Ksp  (26) 

where  I<  is  the  permeability  of  the  porous  media  and  the  constant 
p  takes  the  value  of  3  for  catalyst  layers  and  4.5  for  diffusion  layer 
[23].  In  the  literature,  pc(s)  function  usually  modeled  using  Leverett 
J(s)  function  approach  where  Udell’s  empirical  correlation  used  for 
pc(s)  relationship  [23].  However,  this  correlation  is  not  appropriate 
to  use  in  PEM  fuel  cell  modeling  since  it  is  for  sand/rock  type  porous 
media.  Thus,  in  the  present  model  experimentally  determined  cap¬ 
illary  functions  given  by  Ye  and  Nguyen  [23]  were  used. 

The  source  terms  of  all  governing  equations  of  the  present  model 
are  indicated  in  detail  with  respect  to  PEM  fuel  cell  components 
in  Tables  2-4.  Definition  of  each  of  the  source  terms  given  in 
Tables  2-4  is  defined  in  Table  5. 

The  resistance  to  gas  flow  in  diffusion  and  catalyst  layers  are 
given  by  Eqs.  (27)-(29)  and  gas  permeability  correlation  (1  -s)p 


Table  2 

Source  terms  of  mass  and  momentum  conservation  equations. 


Smass 

Smom-x 

Smom-y 

Smom-z 

Anode 

Mono-polar  plate 

- 

- 

- 

- 

Flow  channel 

- 

- 

- 

- 

Diffusion  layer 

—Sphase 

Sm-x 

Sm-y 

Sm-z 

Catalyst  layer 

Sh2  —  $  phase  ~  S(HssM  h20 

Sm-x 

Sm-y 

Sm-z 

Membrane 

- 

- 

- 

- 

Cathode 

Catalyst  layer 

So2  ^ phase  S^issM  h20 

Sm-x 

Sm-y 

Sm-z 

Diffusion  layer 

—Sphase 

Sm-x 

Sm-y 

Sm-z 

Flow  channel 

- 

- 

- 

- 

Mono-polar  plate 

- 

- 

- 

- 

takes  different  exponents  p  as  shown  in  Table  5  [23].  The  expres¬ 
sions  for  consumption  of  hydrogen  and  oxygen  and  production 
of  water  are  given  in  Eqs.  (30)-(32).  In  the  present  model,  water 
is  assumed  to  be  formed  in  the  ionomer  phase  [27].  Ohmic  heat¬ 
ing  (Eqs.  (33)  and  (34)),  electrochemical  reaction  heat  source  (Eq. 
(35))  and  phase  change  source  expressions  (Eqs.  (36)  and  (37))  are 
defined  as  in  Table  5.  The  source  terms  of  electric  and  proton  poten¬ 
tial  equations  (Eq.  (38))  are  calculated  from  the  volumetric  transfer 
current  (&t).  For  anode  catalyst  layer,  the  volumetric  transfer  current 
was  calculated  by  using  the  Butler-Volmer  equation: 


(1  -  s)Scrpt 


mPt  Tref 


x 


exp 


/  otaFr)a  \ 

V  RT  ) 


exp  (-(1  -afl)^) 


(42) 


where  pa  is  the  anode  over-potential  defined  as  the  difference 
between  the  electric  and  proton  potentials: 


da  —  0e  —  0p  (43) 

For  the  cathode  catalyst  layer  following  expression  was  used  to 
calculate  the  volumetric  transfer  current: 

9?  (Am3)  =  4F(\  —  s)(3/ragg)(gso,  +  e,-agg) -  ^  (44) 

3/(deffaaggk(DRRragg )  +  r agg3/(Do2(ragg  +  5))  Ho2 

where  the  first  term  in  the  denominator  represents  the  resistance 
due  to  diffusion  and  reaction  inside  the  spherical  agglomerates 
whereas  the  second  term  represents  the  resistance  due  to  ionomer 
film  surrounding  the  spherical  agglomerates.  Effect  of  blockage  of 
active  area  by  liquid  water  is  modeled  by  the  term  (1  -s)  in  the 
numerator  of  Eq.  (44)  [28].  The  internal  effectiveness  factor  ( r]eff ) 
for  the  spherical  agglomerates  is  given  as  follows  [29] 

"«-(£)  (sabri)  (45) 


Table  3 

Source  terms  of  individual  species  and  energy  conservation  equations. 


SSp.H2 

$sp.  02 

SSp.H20 

$th 

Anode 

Mono-polar  plate 

- 

- 

- 

S ohm-e 

Flow  channel 

- 

- 

- 

- 

Diffusion  layer 

- 

- 

—Sphase 

Sohm-e  +  $lat 

Catalyst  layer 

Sh2 

- 

phase  —  $dissM  h20 

Sohm-e  +  Sohm-i  +  Srxn  +  $!at  +  ^sd 

Membrane 

- 

- 

- 

Sohm-i 

Cathode 

Catalyst  layer 

- 

So2 

$  phase  $dissM  h20 

Sohm-e  +  Sohm-i  +  Srxn  +  $lat  +  Ssd 

Diffusion  layer 

- 

- 

—Sphase 

Sohm-e  +  Slat 

Flow  channel 

- 

- 

- 

- 

Mono-polar  plate 

- 

- 

- 

Sohm-e 
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Table  4 

Source  terms  of  electric  potential,  proton  potential,  dissolved  water  and  liquid  water 
conservation  equations. 


se 

sp 

Sdw 

5S 

Anode 

Mono-polar  plate 

- 

- 

- 

- 

Flow  channel 

- 

- 

- 

- 

Diffusion  layer 

- 

- 

- 

5 phase 

Catalyst  layer 

—Strn 

+5frn 

Sdiss  Sosm 

5 phase 

Membrane 

- 

- 

Sosm 

- 

Cathode 

Catalyst  layer 

+Strn 

—Strn 

Sdiss  "F  Sosm  +  Sh2  0 

5 phase 

Diffusion  layer 

- 

- 

- 

5 phase 

Flow  channel 

- 

- 

- 

- 

Mono-polar  plate 

- 

- 

- 

- 

where  cp  is  a  dimensionless  group  known  as  Thiele  modulus  [29] 
for  chemical  reactions  and  for  spherical  particles  it  is  defined  as: 


<P  = 


(46) 


where  ragg  is  the  spherical  agglomerate  radius  and  aagg  corre¬ 
sponds  the  active  area  per  unit  agglomerate  volume  (Eq.  (6)).  In 
the  model,  it  was  assumed  that  the  spherical  agglomerates  are 
porous;  therefore,  diffusion  coefficient  of  the  oxygen  in  the  spher¬ 
ical  agglomerates  should  be  modified  accordingly.  In  the  present 
work,  Bruggemann’s  correction  is  employed  in  the  calculation 
of  the  effective  inter-agglomerate  diffusion  coefficient  of  oxygen 
(Do?) 


Do?  =  ^o2 


(47) 


where  Do2  is  the  diffusion  coefficient  of  oxygen  in  the  ionomer. 
The  rate  coefficient  (I<orr)  for  oxygen  reduction  reaction  (ORR)  is 
derived  by  writing  the  rate  of  reaction  of  oxygen  in  first  order 


Table  5 

Definition  of  individual  source  terms. 


Term 

Expression 

Sm-x 

. u 

W-s? 

5  m-y 

.  V 

K'yO-SjP 

Sm-z 

■  ^ g  w 

K'd-sf 

Sh2 

-§Mh2 

5o2 

-§Mo2 

5  h2o 

+  §Mh2o 

5  ohm-e 

*e_x  !e-y  'e-z 

aeff  +  aeff  +  aeff 

s-x  s-y  s-z 

Sohm-i 

Pp/of 

5  rxn 

Slat 

Sphase  he 

Ssd 

Sdiss^H2ohc 

Strn 

ft 

Sphase 

kcS o  (1  -S)P-22!^fM. 
keS0s{Pwv  —  f*sat)PH; 

Sdiss 

kdw.(C 2  -  Qw) 

Sosm 

V-  (n,j3lV^p) 

sat 


2°(l )  TPy 


<  p 


sat 


reaction  rate  form  (-r  =  /<orrQ)2)  and  using  Butler-Volmer  type 
kinetics  [6,14]  for  ORR 


exp  (-^r)-exp  ((1~“c)?r) 

°2 

where  cathode  over-potential  (rjc)  is  given  by 
IJc  =  (pe  —  (pp  —  Vocv 


k<DRR  = 


Jo 


Apr 


ref 


(48) 


(49) 


Open  current  voltage  (Vqcv)  °f  the  PEM  fuel  cell  was  calcu¬ 
lated  using  Eq.  (50)  correlated  from  the  experimental  data  given 
by  Parthasarathy  et  al.  [30]. 


Vqcv  =  0.002534 T  +  0.9251  (50) 

In  Eq.  (44),  represent  the  reference  oxygen  concentration 
at  the  conditions  valid  for  the  reference  exchange  current  density 
( J'jf)  of  oxygen  reduction  reaction.  The  experimental  data  from  the 
series  of  work  published  by  Parthasarathy  et  al.  [30,31]  was  used 
to  correlate  the  reference  exchange  current  density  to  temperature 
and  used  to  evaluate  the  transfer  coefficient.  The  transfer  coefficient 
( ac )  was  calculated  depending  on  the  operating  voltage  of  the  PEM 
fuel  cell  as  follows: 

fl  Vcen  >  0.8  V  1 

c  f  0.495  +  2.3  x  10-3(T  —  300)  Vce//<0.8Vj  1  J 

In  the  diffusion  and  catalyst  layers,  water  may  condense  and 
block  the  pores  in  the  catalyst  layer  if  water  partial  pressure  exceeds 
the  saturation  pressure.  The  difference  between  saturation  and  par¬ 
tial  pressure  of  water  is  taken  as  the  driving  force  for  the  phase 
change.  Thus,  the  phase  change  of  water  is  modeled  using  the 
expression  given  in  Eq.  (39)  where  condensation  (kc)  and  evapo¬ 
ration  (ke)  rates  are  assumed  to  be  constant  [32]. 

In  catalyst  layers,  water  may  also  dissolve  into  the  ionomer 
phase  according  to  water  sorption  equilibrium  of  the  ionomer.  The 


Note 


i  =  DL-*p  =  4.5 

i  =  CL -^p  =  3 

(27) 

i  =  DL^p  =  4.5 

i  —  CL  — >  p  =  3 

(28) 

i  =  DL^p  =  4.5 

i  =  CL ->p  =  3 

(29) 

(30) 

(31) 

(32) 

Reduces  to  i^/cre/  for  isotropic  medium 

(33) 

(34) 

n  =  2  for  anode  CL 

n  =  4  for  cathode  CL 

(35) 

(36) 

(37) 

Eq.  42  and  44 

(38) 

£o  is  the  dry  porosity  of  the  porous  medium. 

(39) 

(40) 

(41) 
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sorption/desorption  of  water  is  modeled  using  Eq.  (40)  where  sorp¬ 
tion/desorption  rate  constant  kdw  is  assumed  to  be  100s_1  [27,33] 
in  order  to  maintain  near  equilibrium  conditions.  The  driving  force 
for  the  sorption/desorption  is  the  difference  between  equilibrium 
concentration  of  dissolved  water  in  the  ionomer  phase  and  its 
actual  value.  In  the  dissolved  water  conservation  equation,  water 
transport  due  to  electro-osmotic  drag  (Eq.  (41 ))  is  treated  explicitly 
as  a  source  term.  In  the  present  model,  electro-osmotic  drag  coeffi¬ 
cient  ( nd )  was  taken  as  a  linear  function  of  water  content  (A.)  of  the 
polymer  electrolyte  [25]: 

-  - 12#) 

and  the  equilibrium  water  content  of  the  ionomer  and  the  mem¬ 
brane  was  modeled  using  the  expression  given  by  Hinatsu  et  al. 

[34] : 

Xeq{a)  =  0.30  +  10.8a  -  16.0a2  +  14.1a3  (53) 

where  a  is  the  relative  humidity  {Ph2o/P^0)  of  the  gas  mixture. 
For  relative  humidity  values  larger  than  unity,  the  equilibrium 
water  content  is  fixed  to  constant  value  of  Xeq(\)  by  neglecting 
the  Schroeder’s  Paradox  since  the  absence  of  it  has  been  illustrated 

[35] . 

2.4.  Numerical  methods 

A  single  channel  PEM  fuel  cell  with  straight  flow  channels 
on  both  anode  and  cathode  sides  was  considered  in  this  work. 
In  the  modeling,  single  domain  approach  [36]  was  used  so  that 
requirement  of  internal  boundary  conditions  was  avoided.  In  this 
approach,  all  of  the  governing  equations  solved  numerically  for  all 
components  of  the  PEM  fuel  cell.  However,  special  treatments  were 
applied  to  diffusion  and  convection  terms  in  order  to  avoid  the 
transport  of  dependent  variable  through  the  PEM  fuel  cell  compo¬ 
nent  on  which  dependent  variable  does  not  physically  exist  [23,36]. 

The  governing  mass,  momentum,  energy,  charge,  dissolved 
water  and  liquid  water  conservation  equations  were  solved  using 
the  finite-volume  technique.  The  commercial  flow  solver  FLUENT® 
was  utilized  with  its  user  defined  function  (UDF)  feature  by  which 
charge,  dissolved  water  and  liquid  water  conservation  equations 
were  implemented.  Diffusion  coefficients  and  source  terms  of  all 
conservation  equations  as  well  as  all  physical  property  defini¬ 
tions,  correlations  and  geometrical  definitions  were  incorporated 
by  using  custom  written  user  defined  C-language  functions. 

The  SIMPLE  algorithm  [37]  was  used  to  handle  the 
pressure-velocity  coupling  in  the  solution  of  continuity  and 
momentum  equations.  Due  to  intrinsic  non-linearity  of  the 
present  model,  the  under-relaxation  and  the  source  term  lin¬ 
earization  techniques  [37]  were  applied.  Detailed  grid  dependence 
study  was  performed  to  maintain  the  solution  time  as  low  as 
possible  while  keeping  the  accuracy  of  the  solution  reasonable. 
For  the  present  study,  the  total  number  of  control  volume  was 
fixed  to  95,040  for  the  half  of  the  single  channel  geometry.  The 
convergence  criteria  of  5  x  10-6  were  used  for  the  total  residual  of 
each  conservation  equation. 

3.  Results  and  discussion 

3.2.  Model  validation 

The  three-dimensional  model  presented  in  preceding  section 
validated  against  the  experimental  data  reported  by  Chang  et  al. 
[38].  The  MEA  in  the  work  of  Chang  et  al.  [38]  has  50  cm2  and  has 
multi-channel  serpentine  flow  field.  In  this  study,  one  of  the  chan¬ 
nels  of  the  PEM  fuel  cell  is  considered.  The  cross  flow  between 
adjacent  channels  and  the  channel  bends  are  neglected.  Thus, 


Table  6 

Key  electrode  parameters  for  baseline  case. 


Anode  electrode  parameters 


Jo 

100  Am-2 

aa 

0.5 

Cre/ 

0.59  x  10-3  kmolm-3 

[39] 

TptC 

0.2 

rPt 

0.75 

[6,22] 

Si 

0.59 

cCL 

Sq 

0.2 

PDl 

£0 

0.74 

[40] 

Cathode  electrode  parameters 

Jo, LCD 

1.3665  x  10-4  Am-2 

[30,31] 

Jo.HCD 

1.8694  x  10-2  Am-2 

[30,31] 

Ea,lcd 

73.3  kcalmol-1 

[30] 

Ea.hcd 

27.6  kcal  mol-1 

[30] 

Cref 

2.281  x  10“3  kmolm-3 

[30,31,39] 

SptC 

0.2 

Baseline  value 

rPt 

0.75 

[6,22] 

ragg 

0.3  |jim 

Baseline  value 

£agg 

0.45 

Baseline  value 

s 

85  nm 

Baseline  value 

oCL 

£0 

0.2 

Baseline  value 

PDL 

£0 

0.74 

[40] 

straight  and  single  channel  fuel  cell  geometry  with  25  cm  length, 
as  stated  by  Chang  et  al.  [38]  is  used  as  model  geometry.  The  fuel 
cell  dimensions,  MEA  parameters  and  component  properties  used 
for  the  validation  is  taken  either  from  the  data  given  in  Chang  et  al. 
[38]  or  taken  from  the  literature  and  the  manufacturer  data  (see 
Tables  6  and  7).  Fig.  2  compares  the  polarization  curve  obtained 
from  the  present  three-dimensional  model  and  the  experimental 
data.  Over  the  whole  range  of  experimental  data,  good  agree¬ 
ment  between  numerical  and  experimental  polarization  curve  is 
obtained.  Small  deviation  between  numerical  and  experimental 
polarization  curves  at  high  current  densities  may  be  attributed 
to  neglected  channel  bends  and  cross-flow  between  the  adjacent 
channels  as  well  as  parameter  uncertainties. 

3.2.  Analysis  of  cathode  catalyst  structure 

After  validating  the  present  model  with  the  experimental  data, 
parametric  analysis  of  cathode  catalyst  layer  structure  were  per¬ 
formed  in  order  to  see  the  changes  in  area  specific  and  platinum 

Table  7 

Properties  and  constants  for  baseline  case. 


Electrical  properties  of  anode  and  cathode  side  fuel  cell  components 


<jdl 

uxy 

17241  Sm-1 

(jDL 

J  z 

1250  Sm-1 

0CL 

300  Sm-1 

aMP 

20000  Sm-1 

kdl 

xy 

1.3  x  10-11  m2 

[41] 

I<f>L 

2.9  x  10-12  m2 

[41] 

I<CL 

3.0  x  10-14  m2 

[23] 

Thermal  properties  of  anode  and  cathode  side  fuel  cell  components 

l<MP 

50  W  m-1  K-1 

kDL 

xy 

11.0  Wm-1  K-1 

[42] 

fcf 

0.88  Wm-1  K-1 

[42] 

1<CL 

0.27  Wm-1  K-1 

[43] 

l<M 

0.1 6  Wm-1  K-1 

[43] 

Membrane  properties 

Pm 

1970  kg  m-3 

EWm 

lOOOkgkmol-1 

Additional  constants 

ahz 

—226  x  103  J  kmol-1  K-1 

[44] 

AHZ 

62.8  x  103  Jkmol-1  K-1 

[44] 

kc 

100  s-1 

[32] 

ke 

9.86  x  10-4  Pa-1  s-1 

[32] 

kdw 

100  s-1 

[27] 
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Current  Density,  Acm  -2 

Fig.  2.  Comparison  of  the  model  prediction  with  the  experimental  data  reported  by 
Chang  et  al.  [38]. 


Table  8 

Geometric  dimensions  of  the  single  channel  PEM  fuel  cell. 


Widtha/height/length  of  flow  channels 

Width3  of  flow  channel  spacing  (shoulder) 

Thickness  of  diffusion  layers 

Thickness  of  anode  catalyst  layer 

Thickness  of  cathode  catalyst  layers  (baseline  value) 
Thickness  of  membrane 

Thickness  of  mono-polar  plates 

0.5  mm/1  mm/50  mm 
0.5  mm 

0.16mm 

0.02  mm 

0.04  mm 

0.04  mm 

1.2  mm 

3  The  value  corresponds  to  half  of  the  fuel  cell  geometry. 

Table  9 

Operating  parameters  for  baseline  case. 

Anode  feed 

Stoichiometry  at  1  Acm-2 

1.5 

Mass  flow  rate 

2.164  x  lO"8  kgs”1 

Temperature 

65  °C 

Relative  humidity 

100% 

Cathode  feed 

Stoichiometry  at  1  Acm-2 

2.0 

Mass  flow  rate 

3.998  x  10“7  kgs-1 

Temperature 

65  °C 

Relative  humidity 

100% 

Anode/cathode  terminal  temperatures 

65  °C 

Anode/cathode  side  back  pressures 

152  kPa 

mass  specific  polarization  and  power  density  curves  with  chang¬ 
ing  cathode  catalyst  parameters.  For  this  purpose  a  baseline  case 
was  prepared  and  its  polarization  curve  was  obtained  by  using  the 
present  model.  The  value  of  each  cathode  catalyst  layer  parameters 
is  doubled  or  halved  while  keeping  the  remaining  parameters  con¬ 
stant.  Polarization  curves  generated  by  the  model  were  compared 
in  order  to  investigate  their  effects  on  the  PEM  fuel  cell  perfor¬ 
mance.  Following  catalyst  layer  parameters  were  investigated  in 
the  present  analysis: 

•  Platinum  mass  ratio  of  catalyst  particles  (rPtc). 

•  Catalyst  layer  thickness  (tcat). 

•  Spherical  agglomerate  radius  (ragg). 

•  Agglomerate  porosity  (£agg). 

•  Ionomer  film  thickness  around  an  agglomerate  (<5). 

•  Catalyst  layer  porosity  (£cat). 


domain  by  dividing  it  through  the  middle  of  the  channel  (Fig.  3). 
Same  flow  channel  and  diffusion  layer  dimensions  are  used  for 
both  side  of  the  fuel  cell.  The  geometrical  dimensions  are  listed 
in  Table  8.  The  operating  conditions  are  as  follows:  fully  humidi¬ 
fied  hydrogen  gas  and  fully  humidified  air  at  65  °C  was  fed  into  the 
anode  and  cathode  sides,  respectively.  The  temperature  of  anode 
and  cathode  terminals  of  the  fuel  cell  was  fixed  to  65  °C.  Anode  and 
cathode  channel  back  pressures  were  set  to  152  kPa.  The  fuel  cell 
electrode  properties,  fuel  cell  parameters  and  the  operating  param¬ 
eters  used  in  the  baseline  simulation  are  listed  in  Tables  6-7  and  9, 
respectively.  The  values  of  catalyst  layer  parameters  used  for  para¬ 
metric  analysis  are  summarized  in  Table  10.  Totally  13  polarization 
curves  were  generated  during  the  simulations,  each  requiring  runs 
for  eight  different  operating  voltages.  In  the  subsequent  sections, 
these  polarization  curves  are  compared  and  discussed  in  detail. 


For  the  baseline  case,  single  channel  PEM  fuel  cell  geometry 
is  considered.  Half  of  the  cell  geometry  is  used  as  computational 


mono¬ 
plate 

flow  eh; 

diffusion 

catalyst 

catalyst  layer 

diffusion  layer 

flow  channel 

mono-polar 
plate 

Fig.  3.  Geometry  of  half  of  the  straight  channel  PEM  fuel  cell  used  in  this  study. 


3.2.2.  Effect  of  platinum  mass  ratio  of  catalyst  particles 

The  platinum  mass  ratio  of  catalyst  particles  determines  the 
platinum  loading  of  catalyst  layer  for  constant  values  of  other 
catalyst  parameters.  For  rPtc  ratios  of  0.1,  0.2  and  0.4,  the  plat¬ 
inum  loadings  of  catalyst  layer  are  approximately  0.46,  1.02 
and  2.61  mgPtcnrr2,  respectively  while  volume  fraction  of  solid, 
ionomer  and  gas  phase  as  well  as  ionomer  loading  (2.2  mg  cm-2) 
are  constant  for  all  cases.  In  addition,  increase  of  rPtc  ratio  also 
decreases  the  specific  active  area  of  catalyst  particles  as  indicated 
in  Eq.  (7).  Thus,  these  two  competitive  effects  determine  the  total 
active  area  per  unit  volume  of  catalyst  layer  that  is  proportional  to 
catalyst  rPtc  ratio  for  fixed  catalyst  layer  thickness.  Hence,  the  fuel 
cell  performance  is  expected  to  be  increased  when  rPtc  increased. 
Fig.  4  shows  the  polarization  and  the  power  curves  obtained  using 
the  present  model  for  platinum  mass  ratio  of  0.1,  0.2  and  0.4.  Con¬ 
sidering  the  change  of  platinum  loading,  polarization  and  power 
curves  (per  unit  area)  of  high  rPtc  ratio  case  shifted  upward  while 


Table  10 

Cathode  catalyst  layer  parameters  used  for  parametric  analysis. 


rptc 

teat  (fxm) 

ragg([xm) 

£agg 

<5  (nm) 

Scat 

Low  value 

0.1 

20 

1 

0.15 

25 

0.1 

Baseline 

0.2 

40 

2 

0.30 

50 

0.2 

High  value 

0.4 

80 

4 

0.60 

100 

0.4 
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Fig.  4.  Effect  of  platinum  mass  ratio  of  catalyst  particles  ( rPtc )  on  the  fuel  cell  perfor¬ 
mance  for  tCat  =  40  |jim,  5  =  50  nm,  ragg  =  2  p.m,  sagg  =  0.3  and  ecat  =  0.2.  (a)  Polarization 
and  power  curves  per  unit  MEA  area  and  (b)  polarization  and  power  curves  per  unit 
platinum  mass. 


Fig.  5.  Effect  of  catalyst  layer  thickness  (tcat)  on  the  fuel  cell  performance  for 
rPtc  =  0.2,  5  =  50nm,  raffi  =  2|xm,  sagg  =  0.3  and  ecat  =  0.2.  (a)  Polarization  and  power 
curves  per  unit  MEA  area  and  (b)  polarization  and  power  curves  per  unit  platinum 
mass. 


that  of  low  rPtc  ratio  shifted  downward  (Fig.  4a).  However,  the  cur¬ 
rent  and  power  output  from  the  fuel  cell  per  unit  mass  of  platinum 
shows  inverse  behavior  so  that  better  performance  obtained  for  low 
rPtc  ratio  (Fig.  4b).  This  implies  that  increase  of  platinum  loading  by 
rPtc  ratio  may  not  be  as  favorable  as  expected  in  terms  of  effective 
utilization  of  the  catalyst  because  of  the  decreasing  specific  surface 
area  of  catalyst  particles.  From  the  current  analysis,  the  increase 
of  rPtc  ratio  results  in  higher  area  specific  current  and  power  den¬ 
sities  for  the  range  of  parameters  studied.  However,  the  effective 
utilization  of  platinum  loading  was  attained  for  low  value  of  rPtc. 

3.2.2.  Effect  of  catalyst  layer  thickness 

Fig.  5  shows  the  polarization  and  power  curves  of  simulations 
obtained  by  altering  the  catalyst  layer  thickness  from  its  base  value 
(40p>m)  to  20  and  80p>m.  The  operating  conditions  as  well  as 
other  catalyst  parameters  are  kept  constant.  The  direct  effect  of 
catalyst  layer  thickness  is  on  platinum  loading  of  catalyst  layer 


of  the  fuel  cell;  the  higher  the  catalyst  layer  thickness  the  larger 
the  platinum  loading  of  catalyst  layer.  From  Fig.  5a,  it  is  seen  that 
for  low  current  densities,  area  specific  current  and  power  output 
per  unit  area  of  cell  increases  with  catalyst  layer  thickness.  This 
is  because  of  lower  oxygen  demand  of  the  catalyst  layer  at  low 
current  densities  resulting  in  excess  amount  of  oxygen  in  the  cat¬ 
alyst  layer.  Therefore,  current  density  per  unit  area  increases  with 
tCat  at  low  current  densities.  However,  at  higher  current  densities 
>0.75  A  cm-2  approximately  the  same  area  specific  current  density 
was  observed  for  40  and  80  p,m.  This  implies  that,  the  determining 
factor  for  power  output  from  the  fuel  cell  starts  to  become  the  lim¬ 
itation  of  mass  transport  inside  the  catalyst  layer.  At  this  point  use 
of  thicker  catalyst  layer  is  ineffective. 

On  the  other  hand,  if  current  and  power  output  of  fuel  cell  com¬ 
pared  on  platinum  mass  basis  (Fig.  5b),  the  thinnest  catalyst  layer 
(20p>m,  ~0.51  mgptcnrr2)  displayed  better  performance  in  terms 
of  effective  platinum  usage  whereas  fuel  cell  performance  in  terms 
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Fig.  6.  Change  of  catalyst  layer  oxygen  concentration  along  y-direction,  at  the  mid¬ 
dle  of  catalyst  layer  inx-direction  and  10  |jim  away  from  the  membrane  at  0.4  V  cell 
potential. 

of  per  unit  area  is  lower  than  its  base  case  (tcat  =  40  |jim).  One  of  the 
major  reasons  for  that  is  the  resistance  to  oxygen  transport  in  cat¬ 
alyst  layer.  In  the  thinnest  cathode  catalyst  layer,  average  oxygen 
concentration  in  the  cathode  catalyst  layer  is  expected  to  be  higher. 
This  is  due  to  the  shorter  average  transport  length  of  oxygen  from 
diffusion  layer  to  catalyst  active  sites.  This  effect  is  clearly  seen  in 
Fig.  6,  where  the  oxygen  concentration  along  the  channel  direction 
(y-direction)  is  plotted  for  the  same  distance  away  the  membrane. 
According  to  Fig.  6,  the  oxygen  concentration  along  the  channel  is 
always  larger  for  20  |jim  thick  catalyst  layer.  Another  reason  might 
be  the  decrease  of  cathode  over  potential  because  of  the  thin  cat¬ 
alyst  layer  that  means  lower  total  resistance  to  proton  transport. 
Consequently,  for  the  range  of  interest,  even  though  thicker  cata¬ 
lyst  layer  resulted  in  higher  current  and  power  output  per  unit  area, 
effective  utilization  of  catalyst  layer  might  be  obtained  by  using 
thinner  catalyst  layer. 

3.2.3.  Effect  of  spherical  agglomerate  radius 

In  the  spherical  agglomerate  model,  as  the  agglomerate  radius 
increases  at  constant  catalyst  thickness  and  porosity;  solid  vol¬ 
ume  fraction  increases  and  ionomer  volume  fraction  decreases  as 
indicated  by  Eq.  (3).  Thus,  platinum  loading  of  catalyst  layer  and 
agglomerate  radius  are  proportional.  If  the  change  of  platinum 
loading  with  the  agglomerate  radius  is  considered,  the  performance 
of  the  cell  should  increase;  however,  in  the  present  model,  inter¬ 
particle  effectiveness  factor  of  spherical  agglomerates  are  inversely 
proportional  with  ragg  as  pointed  in  Thiele  modulus  (Eq.  (46)).  In 
addition,  the  increase  of  agglomerate  radius  also  enhances  electron 
transport  by  positive  contribution  to  effective  electrical  conduc¬ 
tivity  of  catalyst  layer  (Eq.  (21)).  At  the  same  time,  reduction  in 
ionomer  volume  fraction  adversely  affects  effective  proton  con¬ 
ductivity  of  catalyst  layer  so  that  effective  proton  conductivity  of 
catalyst  layer  decreases  (Eq.  (22)). 

Polarization  and  power  curves  for  three  different  values  of 
agglomerate  radius  are  shown  in  Fig.  7.  It  is  observed  that  current 
and  power  output  per  unit  area  and  per  unit  platinum  mass  both 
increase  with  decreasing  agglomerate  radius.  It  is  also  observed 
that  cell  performance  is  influenced  considerably  by  agglomerate 
inter-particle  effectiveness  factor.  The  effect  can  also  be  inferred 


Fig.  7.  Effect  of  spherical  agglomerate  radius  (ragg)  on  the  fuel  cell  performance  for 
rPtc  =  0.2,  <S  =  50nm,  tcat  =  40  |xm,  eagg  =  0.3  and  scat  =  0.2.  (a)  Polarization  and  power 
curves  per  unit  MEA  area  and  (b)  polarization  and  power  curves  per  unit  platinum 
mass. 

from  volume  average  values  of  agglomerate  inter-particle  effec¬ 
tiveness  factor,  effective  electrical  and  proton  conductivity  of 
cathode  catalyst  layer  and  average  current  density  of  the  fuel 
cell.  In  Fig.  8,  normalized  values  of  these  quantities  are  plotted 
against  the  agglomerate  radius  at  0.4  V  cell  potential.  Fig.  8  indi¬ 
cates  that,  the  effective  electrical  conductivity  increases  with  ragg 
and  effective  proton  conductivity  and  inter-particle  effectiveness 
factor  have  reverse  trend.  Besides,  inter-particle  effectiveness  fac¬ 
tor  is  found  to  be  more  sensitive  to  ragg  than  the  other  parameters. 
Therefore,  current  analysis  suggests  that  better  cell  performance 
might  be  obtained  by  increasing  both  the  effectiveness  factor  and 
the  effective  proton  conductivity  by  reducing  the  agglomerate 
radius,  ragg. 

3.2.4.  Effect  of  agglomerate  porosity 

In  the  present  model,  it  was  assumed  that  voids  in  the  spheri¬ 
cal  agglomerates  filled  only  with  ionomer.  The  oxygen  reaches  the 
active  reaction  sites  by  diffusing  through  the  ionomer.  Thus,  larger 
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Fig.  8.  Effect  of  agglomerate  radius  ( ragg )  on  normalized  volume  average  values  of 
inter-particle  effectiveness  factor  ( rjeff ),  effective  electric  and  proton  conductivity  of 
catalyst  layer  op and  average  current  density  (Javg )  at  0.4  V  cell  potential  for 
fixed  values  of  rPtc  =  0.2,  <5  =  50  nm,  tcat  =  40  |xm,  sagg  =  0.3  and  scat  =  0.2. 

agglomerate  porosity  (sagg)  promotes  the  oxygen  diffusion  inside 
the  agglomerate  by  increasing  the  effective  diffusion  coefficient 
given  by  Eq.  (47). 

Fig.  9  shows  the  polarization  and  power  curves  for  £agg  varying 
from  0.15  to  0.6.  The  polarization  curves  in  Fig.  9a  show  that  the 
performance  of  the  PEM  fuel  cell  increases  proportionally  with 
£agg\  however,  at  low  current  densities  <~0.15  A  cm-2  the  perfor¬ 
mance  of  three  cases  are  nearly  the  same  and  the  lowest  obtained 
when  £agg  is  small.  At  high  current  densities,  better  performance 
observed  for  high  value  of  sagg.  The  difference  between  the  three 
cases  might  be  the  effect  of  proton  conductivity  of  cathode  catalyst 
layer  that  increases  with  ionomer  fraction  of  the  catalyst  layer. 
The  values  of  ionomer  volume  fractions  of  three  cases  conform  the 
effect  (ionomer  volume  fractions  are  0.169,  0.280  and  0.503  for 
£agg  of  0.15,  0.3  and  0.6,  respectively).  Therefore,  it  can  be  said  that 
the  PEM  fuel  cell  performance  is  a  strong  function  of  the  ionomer 
volume  fraction  and  so  the  effective  proton  conductivity  of  the 
catalyst  layer. 

On  the  other  hand,  if  the  polarization  curves  were  plotted  in 
terms  of  platinum  loading,  the  relationship  between  sagg  and  the 
cell  performance  become  clearer.  As  shown  in  Fig.  9b,  the  cell 
performance  improves  with  sagg  for  all  current  densities.  This  is 
mainly  due  to  the  fact  that  the  platinum  loading  of  the  catalyst 
layer  decreases  with  increasing  sagg  (Eqs.  (2)-(4)).  The  reduction 
in  proton  transport  resistance  and  platinum  loading  both  increases 
the  mass  specific  current  and  power  output  of  the  cell.  In  addi¬ 
tion,  increase  of  sagg  is  also  beneficial  for  inter-particle  diffusion 
of  oxygen,  since  effective  diffusion  coefficient  of  oxygen  inside 
the  agglomerates  gets  higher  with  sagg ,  thereby  increasing  the 
inter-particle  effectiveness  factor  defined  in  Eq.  (45)  by  means  of 
reducing  the  Thiele  modulus  Eq.  (46).  Therefore,  higher  current 
and  power  density  (per  unit  area  and  per  platinum  loading)  can 
be  obtained  by  large  values  of  sagg. 

3.2.5.  Effect  of  ionomer  film  thickness  around  an  agglomerate 

In  the  present  agglomerate  model,  each  one  of  the  agglom¬ 
erate  is  assumed  to  be  covered  with  a  constant  finite  thickness 


Fig.  9.  Effect  of  agglomerate  porosity  ( eagg )  on  the  fuel  cell  performance  for  rPtc  =  0.2, 
<5  =  50nm,  tcat  =  40  p.m.  ragg  =  2  |jim  and  ecat  =  0.2.  (a)  Polarization  and  power  curves 
per  unit  MEA  area  and  (b)  polarization  and  power  curves  per  unit  platinum  mass. 

ionomer  film.  This  ionomer  film  serves  as  a  connection  between 
agglomerates  for  the  transport  of  protons  through  the  catalyst 
layer.  Moreover,  oxygen  should  diffuse  through  the  ionomer  film 
in  order  to  reach  the  active  surface  of  catalyst  particles.  Thus,  as 
ionomer  film  thickness  (5)  increases  it  will  effect  proton  conduc¬ 
tivity  positively  through  Eq.  (3 );  however,  at  the  same  time  it  causes 
additional  mass  transfer  resistance  (Eq.  (44))  for  transport  of  oxy¬ 
gen  from  catalyst  pores  to  agglomerate  surface. 

In  Fig.  10,  the  polarization  and  power  curves  for  three  different 
values  of  8  (25,  50  and  lOOnm)  are  plotted  for  constant  values  of 
remaining  catalyst  parameters.  Fig.  10a  indicates  that  at  low  cur¬ 
rent  densities  performance  of  the  fuel  cell  is  similar.  Conversely, 
at  high  current  densities  (>0.7  A  cm-2)  the  fuel  cell  with  higher  <5 
(lOOnm)  suffers  from  mass  transport  losses  caused  by  high  value 
of  ionomer  film  thickness.  This  effect  is  illustrated  graphically  in 
Fig.  1 1  in  which  percentage  of  the  transport  resistance  caused  by 
the  ionomer  film  are  given  as  a  function  of  operating  current  den¬ 
sity  for  different  film  thicknesses.  As  shown  in  Fig.  1 1 ,  for  all  values 
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Fig.  10.  Effect  of  ionomer  film  thickness  (5)  on  the  fuel  cell  performance  for  rPtc  =  0.2, 
teat  =  40  |jim,  ragg  =  2  |jim,  eagg  =  0.3  and  ecat  =  0.2.  (a)  Polarization  and  power  curves  per 
unit  MEA  area  and  (b)  polarization  and  power  curves  per  unit  platinum  mass. 


of  film  thicknesses,  its  contribution  to  total  mass  transfer  resis¬ 
tance  is  relatively  small  at  current  densities  lower  than  0.5  Acnrr2 
which  results  in  indistinguishable  fuel  cell  performance  in  this 
region.  However,  at  higher  current  densities  (>0.5  A  cm-2)  trans¬ 
port  resistance  due  to  polymer  film  starts  to  increase  quickly  for 
larger  film  thickness.  It  reaches  nearly  the  half  of  the  total  transport 
resistance  for  the  highest  film  thickness  thereby  reducing  the  per¬ 
formance  of  the  cell.  Meanwhile,  for  film  thicknesses  <50  nm,  the 
resistance  of  the  ionomer  film  is  negligible  even  at  higher  current 
densities. 

Fig.  10b  shows  that  the  fuel  cell  performance  as  current  and 
power  output  per  platinum  loading  increases  with  ionomer  film 
thickness  starting  from  low  current  densities.  This  performance 
improvement  with  <5  is  due  to  increase  of  total  volume  fraction  poly¬ 
mer  phase  (Eq.  (3))  and  decrease  of  platinum  loading  of  cathode 
catalyst  layer  when  other  catalyst  parameters  are  constant.  These 
results  imply  that  increase  of  ionomer  film  thickness  (<$)  might 


Fig.  11.  Variation  of  average  film  resistance  percentage  in  the  cathode  catalyst  layer 
with  polymer  film  thickness  at  various  operating  current  densities  for  rPtc  =  0.2, 
teat  =  40  |xm,  ragg  =  2  |xm,  sagg  =  0.3  and  scat  =  0.2. 


cause  increased  mass  transport  resistances  at  higher  current  den¬ 
sities  per  unit  area;  however,  at  the  same  time,  both  the  increase  of 
total  polymer  phase  volume  fraction  and  the  decrease  of  platinum 
loading  might  significantly  enhance  catalyst  utilization  as  indicated 
by  the  polarization  curves  plotted  for  per  unit  mass  of  platinum 
loading. 

3.2.6.  Effect  of  catalyst  layer  porosity 

Fig.  12  shows  the  variation  of  the  performance  of  the  fuel  cell 
with  cathode  catalyst  layer  porosities  of  0.1,  0.2  and  0.4.  The  key 
effect  of  porosity  is  on  mass  transport  of  oxygen  through  cata¬ 
lyst  layer  pores  that  is  controlled  by  effective  diffusion  coefficient 
defined  in  Eq.  (20).  Thus,  as  porosity  increases  the  diffusion  of 
oxygen  enhances,  vice  versa.  On  the  other  hand,  the  ionomer  and 
platinum  loading  as  well  as  solid  and  ionomer  volume  fractions 
of  catalyst  layer  decreases  as  porosity  of  catalyst  layer  increases 
for  constant  value  of  the  catalyst  layer  thickness.  This  change 
has  an  impact  on  both  electron  and  proton  transport  in  the  cat¬ 
alyst  layer  through  effective  conductivity  values  determined  by 
solid  and  ionomer  void  fractions  of  the  catalyst  layer  (Eqs.  (21) 
and  (22)). 

According  to  Fig.  12a,  once  the  porosity  of  the  catalyst  layer 
increases  the  current  and  power  densities  per  unit  area  of  the  fuel 
cell  deteriorate.  This  implies  that  the  effect  of  effective  electron  and 
proton  conductivities  is  more  significant  than  the  effect  of  oxy¬ 
gen  diffusion  coefficient  resulting  in  better  performance  for  low 
values  of  scat.  Although  the  performance  of  the  fuel  cell  per  unit 
area  decreases  with  scat,  opposite  trend  observed  for  current  and 
power  densities  per  unit  mass.  In  this  case,  the  improvement  of 
oxygen  transport  in  the  catalyst  layers  with  scat  counterbalance 
the  reduction  of  platinum  loading  with  scat  so  that  the  perfor¬ 
mance  of  the  cell  increases  (Fig.  12b).  Therefore,  for  the  present 
values  of  operating  and  catalyst  parameters,  larger  catalyst  layer 
porosity  resulted  in  effective  utilization  of  the  catalyst  loaded,  thus 
larger  mass  specific  power  densities,  however,  at  the  same  time 
lower  current  and  power  output  per  unit  area  obtained  when  scat  is 
high. 
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Fig.  12.  Effect  of  the  catalyst  layer  porosity  (ecat)  on  the  fuel  cell  performance  for 
rPtc  =  0.2,  tcat  =  40  |jim,  5  =  50  nm,  ragg  =  2  pim  and  eagg  =  0.3.  (a)  Polarization  and  power 
curves  per  unit  MEA  area  and  (b)  polarization  and  power  curves  per  unit  platinum 
mass. 

4.  Conclusions 

In  this  paper,  a  three-dimensional,  non-isothermal,  two-phase 
PEM  fuel  cell  computational  fluid  dynamics  model  was  presented 
and  validated  against  the  published  experimental  data.  Ionomer 
film  covered  spherical  agglomerate  model  was  assumed  for  cath¬ 
ode  catalyst  layer.  The  effect  of  catalyst  layer  parameters  such 
as  platinum  mass  ratio  of  catalyst  particles  ( rPtc ),  catalyst  layer 
thickness  ( tcat ),  spherical  agglomerate  radius  (ragg),  ionomer  vol¬ 
ume  fraction  inside  an  agglomerate  (eagg),  ionomer  film  thickness 
covering  the  agglomerates  (5)  and  catalyst  layer  porosity  ( scat )  on 
the  PEM  fuel  cell  performance  were  studied  on  three-dimensional 
geometry  and  for  non-isothermal  PEM  fuel  cell  operation.  Accord¬ 
ing  to  the  numerical  results,  the  effect  of  these  catalyst  layer 
parameters  on  diffusion  coefficients,  electrical  and  proton  con¬ 


ductivities  as  well  as  on  effectiveness  factor  determines  the  area 
specific  and  platinum  mass  specific  current  and  power  densities 
of  the  PEM  fuel  cell  separately.  Thus,  in  the  optimal  design  of  a 
PEM  fuel  cell  catalyst  layer,  one  should  first  choose  which  type  of 
power  output  (area  specific  or  mass  specific)  is  to  be  optimized. 
Then,  depending  on  the  requirements  specified,  optimal  perfor¬ 
mance  can  be  achieved  by  appropriately  adjusting  the  parameters 
of  the  cathode  catalyst  layer. 
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